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The altitude lines of a triangle are concurrent. The bisectors of 
the angles of a triangle are concurrent. Ceva’s theorem says that, all 
the three lines in a triangle which contain a vertex and a point on its 
opposite side are concurrent if and only if no two among them are 
parallel and the product of the three ratios of division of the sides 
made in one direction around the circumference of the triangle is one. 
In Gergonne’s theorem, the three lines of a triangle which are made 
by the vertices and the points of tangency of the incircle on the side 
opposite to them are concurrent. The intersection between the three 
lines tangent to the circumcircle of a triangle and the sidelines opposite 
to them are collinear. 


A point is an extreme point of a plane convex set s unless it lies 
in a triangle which has vertices in s but is no vertex of the triangle. A 
ray from inside a bounded convex figure intersects the boundary of the 
latter at exactly one point. Consecutive vertices of a convex polygon 
exist in sorted angular order about any interior point. A subfacet of a 
simple polytope is shared by two and only two facets. Two facets share 
a subfacet if and only if the latter is determined by d—1 vertices in 
their set; these two facets and the subfacet are called adjacent. A line 
segment defined by two points is an edge of the convex hull if and only 
if all other points of the set lie on, or to one side of it. 


The diameter of a convex figure is the largest distance between 
parallel lines of support. The diameter of its convex hull determines the 
diameter of a set. Every vertex of the Voronoi graph is the intersection 
of three of its edges. Every nearest neighbour of a Voronoi polygon 
defines an edge. 


VT and the triangulation of its nuclei are dual to each other. A 
Voronoi graph on n points has at most 2n —5 vertices and 3n — 6 edges. 
The convex hull of a Voronoi graph on n can be found in linear time. 


Modern programming philosophy puts much emphasise on modu- 
larity of a programme and on information hiding of modules. Though 
undoubtedly information hiding can be good for the finished products, 
during the course of development it sometimes works against yourself 
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when you try to pinpoint an error in order to debug. Some modularisa- 
tions are more about hierarchies than simplicity. Whenever this is the 
case, it is necessary to unconventionally seek a simpler path. 


The explicit equation in 2-d is y = mx+c or ax+by +e = 0. 
Imposing the constraint a? + b? = 1, that is multiplying all the terms by 
(a? + *)-'/?, puts the equation into the canonical or normalised form. 
This makes a = cosa, b = cos 8 and c = —r, where a and bare directional 
cosines, i.e. the cosines of the angles that the normal line makes with 
the x and y axes respectively. Examples of possible conventions are to 
have a normal line point towards outside of the region, to have the line 
direction always to the right of the normal vector, or to keep c positive 
always. 


A parametric form of line equation in 2-d is by introducing a third 
variable ¢ and write the equations as = wa) + ft and y = yo + gt, 
where (20, yo) is the point on the line corresponding to ¢t = 0. A line 
through a point p which makes angles a and 8 with the x and the 
y axes respectively has the parametric equations x = z, + tcosa and 
y = yptcosf. One convention is to vary t from 0 to 1 over a line 
segment, another is to normalise it by multiplying its coefficient by 
(Cig +g?) oe, 


An implicit line equation az + by + c = 0 can be turned into a 
parametric form as 2 = —ac/(a? +b7)'/? + bt and y = —be/(a? +b?)!/? — at. 
And a parametric line described by = a) + ft and y = yo + gt is 
converted into the implicit form as + fy + (xog — yof) =0 


The implicit plane equation in three dimensions is az + by + cz + 
d = 0. The parameters can be found by using Cramer’s rule, a = 
det(1,y;, 2), 6 = det(a;,1,2;), ¢ = det(x;,y;,1) and d = det(2xj, yi, 2), 
which gives a = y1232 + y2213 + y3z21, b = 2 %32 + 22%13 + 23021, € = 
v1y32+ 22413 + x3y21 and d = x1 (y223—y322) + @2(ys21—y123) +23 (y122—y221). 


A normalised form has a constraint a* + b* +c? = 1. This amounts 
to multiplying its implicit equation by (a? +0?+¢?)~!/ to get ax + By + 
yz+6=0. Here a, 8 and ¥ are cosines of the angles which the normal 
to the plane makes with the coordinate axes. The distance between 
two parallel normalised planes is 62 — 6;. A normalised implicit plane 
equation can be used to represent a planar half-space by multiplying 
every tems by —1 and then assign a convention that the vector formed 
by the direction cosines always points towards the outside or the inside 
of the region. 
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The distance from a point to a plane, if the plane is ax+by+cz+d = 
0 and the point is (2p, yp, Zp), is 


a (axp + bYyp + CZp +d)? le 


(a? + b? + c?) 


The intersection of two planes, from the planes ai2 + by +12 + d; = 0 
and at+hytazt+td =0,istx=aot+ ft, y=yotgtand z= 2+ht 
where f = det(b;,¢;), g = det(c;, a;) and h = det(a;, bj), i = 1 and 2. 


The intersection of three planes is found by Algorithm 1, Here the 
minor matrices a’) is 692, 6%£, or 585 as the case may be. 


Algorithm 1 Intersection among three planes. 
ae DL (a,b,0)(-1) aia"; 
if |A| <« then 
at least two of the planes are parallel; 
else 
2 (bide — id's — e188) /d; 
y & (d1d35 — a1 d3§ — e138) /A; 
2 (b,694 + a, 648 — d 58°) /A; 
endif = 
The intersection between a line and the plane az + by +cz +d =0 


is (1 + 2124, y1 + y124, 21 + 212a), where a = —(ax1 + by: +21 + d)/(ax12 + 
by12 + cz12), 12 = 22 — x1 and similarly for yy and 212. 


The area of a circle is mr? and that of its segment is @r7/2. A 
segment is its pie cut reaching its centre while a sector is a plane slice 
through the sphere. The area of a sector is this area subtracted by that 
of a triangle, or r7(6 — sin @)/2. The centre of gravity or the centroid lies 
on the bisector of the central angle with the distance of 4r sin(@/2)/30 
for a sector and 4r sin?(6/2)/3(@ — sin@) for a segment. 
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The volume of a pyramid is Ah/3, where A is the area of base and 
h is the height of the pyramid. The volume of a sphere is 4ar°/3, and 
the distance from its centroid to the sphere centre is 


sin*(6/2) 
Cr/4) | G3 cos(6/2) + cos5(8/2) 


That of a sector of a sphere is (mr? /3)(2—3 cos(6/2) +cos°(/2)) The vol- 
ume of a tetrahedron is V = (1/6) det(x12, #13, £14; Y125 135 Y14) 212) 213) 214), 
where aj; = «;—«4,0r V = (1/6) det(a;, y;, 2,1). The former is limited to 
the case of three dimensions, and is in fact V = (1/6)(a x b)-c, where a, 
b and ¢ are respectively the lines from O to A, B and C in a tetrahedron 
OABC. 


Generalising the latter to higher dimensions, we have the volume 
of a d-dimensional simplex V = (1/d!) det(#;;,1), where 24; is now (xj), 
1<i<(d+1) and1<j<d. Here the information found in existing 
literature seems to be wrong, some lists the multiplying factor as 1/d, 
some simply uses 1/6 throughout all (d > 3)! Tiyapan (2004) mentions 
in more detail how I arrived at the value used here. 

Some of of the algorithms found in literature are the following. 
The algorithm to find whether a point is inside a polygon, Algorithm 2. 
The arbitrary line / here, which is taken for simplicity to be horizontal, 
passes through z. 


Algorithm 2 Point inside a polygon. 


r <Q; 
for i= 1 ton do 
if edge i and / not parallel then 


if i intersects | to the left of z at any point except its lower 
extreme then 


rer+; 
endif 


endif 
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if r odd then 
z is internal to p; 
else 
z is external; 
endif 
endfor cj 
To find the inclusion in a convex polygon, q € p being a known 
fixed point within the polygon, find the wedge in which z lies by doing 


a binary search and test whether Z(zgp;41) is a right turn- while Z(zqp;) 
a left turn angle. If Z(pipi41z) is a left turn angle, then z is inside p. 


The Euclidean minimum spanning tree may be obtained by Algo- 
rithm 3. 


Algorithm 3 Euclidean minimum spanning tree. 


£0; 
for i from 1 to n do 
8(pi) < 0; 
fo pir 
endfor 
while f contains more than one number do 
t¢ f; 
if (s(t)=j) then 


clean up; 
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jegr+h 
endif 
(u,v) + shortest unselected edge incident on t, u € ¢; 
t' ¢ tree in f containing v; 
t’ + merge (t,t')]; 
delete (¢') from f; 
s(t”) + min(s(¢), s(@’)) +1; 
fet’ 
endwhile o 
To rotate v = (v1,v2)' to v' = (vi,v5)1, use v' = Av where A = 


[cos 0, — sin 6; sin 8, cos 6] is the transformation matrix and 6 is the anti- 
clockwise angle of the rotation. 


To rotate a general line in three dimensions by @ around an arbi- 
trary axis, the transformation matrix becomes A = T-!RT, where 
R(6) = Rz'(a)Ry'(8)Rz(0)Ry(8)Re(a)T 
= Ry(—a) Ry(—=B)Rz(8)Ry(B)Re (a) 
The translation T translates one point of the line to the origin and R, 
rotates around the z-axis by a which puts u on to the xz axis, Ry around 
y-axis by 8 and puts wu’ on to the z axis, and R, around z-axis. Since 
the vector v is (x12, y12, 212), we have the a, b, c and the unit vector 
u = (a,b,c)/|v| = (212, y12, 212)/|v|.. Then it follows that cosa = c/d, 
sina = b/d, cos = d and sinB = —a. All of these can perhaps be 
summarised as a linear procedure in Algorithm 4. 


Algorithm 4 Rotation in three dimensions 


v & (212, 912, 212); 


m€ lvl; a a12/m; 
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be yio/m; 
Ce 242/m; 
de (P+); 
Ter + cfd; 
Tay < b/d; 
Tye < d; 
Tyy — —G; 
T < [8(3), —(#1, 1, 21)7;0(3)", 1); 
Ry € [1,0(3)7;0, rae, =T ey; 0:0 jug; 7eas0F0G)", 1]; 
Ry © [Pyaj 0; —Pyy; 02 0,150,003 Sr yy; 0 yx, 0708)", 1); 
R, < [[cos 6, — sin 6; sin 0, cos 6], 0(2, 2); 0(2, 2), $(2)]; 


RT "B5'R,RyRe; v' = Rv. o 


Also in three dimensions, the rotation around the z-axis is 


R, = [{r = 1, m23 = [cos 6, — sin 0; sin 6, cos 6]}] , 


around y-axis is 


Ry = [{r22 = 1,113 = [cos 9, sin 6; — sin 6, cos 6] }| 


and around z-axis is 


R, = [{r33 = 1, m12 = [cos 6, — sin 6; sin 9, cos 6]}] . 


The minor containing parts of the i*® and j*® rows and columns is mj. 


The right hand coordinate system is where a 90° rotation around 
the x-, y- and z-axis bring respectively the y- to z-, z- to g- and g- to y- 
axis. Scaling and translating a vector v in three dimensions amounts to 
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calculating [v’; w] = A[v; w], where A is respectively [$(3)s,0(3);0(3)7, 1] 
and [$(3), Av;0(3)?, 1], Av = (Az, Ay, Az)’. 


A quaternion can be described as a pair (s,v) of a scalar s and a 
vector v = (a,b,c). The rotation by @ around an axis in the direction of 
a unit vector u is then the quaternion (cos 6/2,usin9/2). Let q = (s,v). 
Then q7! = (s,—v). The multiplication of quaternions is giq = (81,01) - 
(82,02) = (8182 — U1 - U2, 81U2 + $201 + V1 X U2), Where the cross product is 
described in minors as v1 X v2 = [64%%; —67*; 674]. 


Let g represent a rotation. Then a vector p is rotated to p' by P' = 
qPq"', where P and P’ are respectively (0,p) and (0,p’). In simplified 
words, this means p' = s*p+ (p-v)v + 2s(u x p) +u x (v x p). Then 
we have the transformation matrix for the general rotation around wu in 
three dimensions, 


(1—2b°-—2c*)  (2ab — 2sc) (2ac + 2sb) 
Ri(6) = (2ab+2sc) (1 —2a* — 2c’) (2be — 28a) 
(2ac — 2sb) (2be+2sa) (1 —2a* — 267) 


where s = cos@/2 and v = (a,b,c) = usin@/2. Furthermore, if q is a 
rotation by 6; around v; and likewise q. by 4) around v2, then g3 = mm 
is a rotation by 63 = 2cos~!s3 around v3 such that sin 63 > 0. 


Quaternion is an extension of complex number to higher dimen- 
sions where there are three imaginary parts instead of one. It is de- 
fined as gq = s+ia+jb+ kc, where a, b, c and s are real numbers, 
? =f =k? =—-1and ij = —ji=k. 


Let a plane be described by (uv — p)-n = 0, where p is a point on-, 
and n a perpendicular to the plane. This means that, for all points v 
lying in the plane (p,n), (v—p)-n = 0. If the plane (p,n) is transformed 
into (Ap,m), then Av lies in (Ap,m) and consequently A'm = n or 
m = (A‘)“1n. 


A vector normal to the surface remains normal if it is transformed 
by (A‘)7!. When a transformation matrix A has the property A = 
Cems for example rotation, all surface normals remain normal. But 
in general this equality does not hold, so we have for instance the 
non-uniform scaling where these normals cease to be normal. 


To test for the intersection between a ray and a triangle using 
Pliicker’s coordinates is described in Algorithm 5. 


Algorithm 5 Ray and triangle intersection. 
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find Pliicker’s coordinates for vertices and the ray; 
test ray against each of the edges; 


if ray hits an edge, passes all of them clockwise or all of them 
counter-clockwise then 


ray intersects the triangle; 
else 
ray and triangle intersect not; 


endif a 


If these vertices are v1, v2 and v3, and the ray is r = 112 = 2 — 11, 
where r; and rz are any two points on the ray, then Pliicker’s coor- 
dinates for the vector v2 are (u,v) = (v2 — v1,U2 x v1), and similarly 
for v23, U31 and 112. For the test between the ray and each edge, find 
C= U,- Uj; +, - uz. Then the ray r counter-clockwise passes the edge i, 
hits it or clockwise passes it respectively asc <0,c=0orc>0. 


A 3-d line can be represented by the six numbers that come with 
coordinates of two distinct points, or by the eight numbers that come 
with the coordinates of two distinct planes. Plticker’s coordinates, how- 
ever, provides a mean which suits geometrical computation better than 
both of these. It redefines the coordinates as u = p—q andv=px4q, 
where p and gq are two points on a line, neither quantity of which de- 
pends on p or g. 


Let a tetrahedron has its vertices at a;, i = 1 to 4. Then the centre 
of its circumsphere is at a; + 6 and its corresponding radius r = |6|, 
where 


6 = [lar2|?(a13 x ara) + Jars|?(ara x @i2) + lara|?(a12 x @i3)] /2|aj2; af3; at4|- 


A triangle Apipop3, where pj; = (a;,y;) and i = 1 to 3, has an 
area A = |z,y,1(3)|, which is positive if and only if Apipop3 forms a 
counter clockwise cycle or /pip2p3 is left-turned. Or equivalently the 
area is A = det(%12, 213; y12, y13), which is positive if the points are in 
anti-clockwise order of the indices and negative otherwise. Or the area 
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is A = [s(s — d,)(s — d))(s d;))'?, where the the semi perimeter s is 
s = )0,d;/2, where d;, i = 1 to 3, are the lengths of the three sides of 
the triangle. 


For a convex polygon, the area can be found by adding together 
the n triangles formed by any two adjacent vertices and one fixed 
point within the polygon. Here n is the number of vertices it con- 
tains. Or it can be found by adding the (n — 2) triangles formed by 
one fixed vertex and any two adjacent vertices of those remaining. An- 
other way of finding the area is A = (1/2) Pa (xiYit1 — YiXi41), Where 
(x;,y;) are vertices. Rearranging to make it faster and more accurate, 
A = (1/2) ar (a4 + 2i+1)(Yir1 — ys)). If the dimension of the poly- 
gon is higher than two, A = (1/2) |N- ™)' (wi x vist) where N is a 
unit vector normal to the plane. The area of a polygon can also be 
computed, without loosing generality, by subtracting the area under its 
lower edges by that of its upper edges. 


A vector normal to a plane is simply n = v2 x v3. 


The ordering of vertices on a face of a polygon is done for the 
purpose of drawing it or for finding vertice pairs which form edges. 
This can be done in two ways. One is to get inside the polygon and look 
at all the vertices around comparing their angles relative to one another. 
Another one is to look at the polygon from a distance and compare their 
angles as before, as well as their distance from the viewing point. The 
angle is 6 = arccos[(v1 - v2)/(|vi||v2|)], where v1 and v2 are vectors to 
vertices from the distant point, and @ the angle between them. 


To find the convex hull in three dimensions, one may use Algo- 
rithm 6. 


Algorithm 6 Convex hull in three dimensions. 
sort s by 2; such that x;(p;) < a;(p;) if and only if i < J; 
if |s < k| then 
construct c(s); 


else 


$1 {p1, --+>Pl[n/2| te 
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82 + {PIn/2|>--->Pn}i 
pi = ($1); 
Pz = ($2); 
p< merge p; and py; 


endif a 


Here c(s) is the convex hull of s. Two convex hulls are merged 
with each other by first constructing a cylindrical triangulation T which 
supports p; and pz along two circuits e; and e2 respectively, then remove 
from both p; and p) the portions which have been obscured by T. 


The following Jarvis’s march algorithm, Algorithm 7, finds a con- 
vex hull in two dimensions. 


Algorithm 7 Convex hull in two dimensions 


pi + the lowest point in s; 
q + the highest point in s; 
while next point # qi do 


find p; € s, 1 = 2,3,..., with an increasing order of the polar 
angles 


with respect to pi; 
endwhile 
while next point 4 p; do 


find q; € s, i = 2,3,..., with an increasing order of the polar 
angles 


with respect to q; and the negative x axis; 
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endwhile a 

The polar angle is an angle with respect to the positive z-axis. The 
lowest and the highest points are on the convex hull. 

Algorithm 8 is the quick hull algorithm. 


Algorithm 8 Quick hull algorithm. 


1 (x0, yo); 
r € (x0, yo + €); 
if s = {I,r} then 
return (1,1); 
else 
find k € s that gives max Aagir or (max Aagir and max Lkir); 
s! < p€s, such that pis on the left of Th; 
s* + q €s, such that gq is on the left of irs 
{h} € (s!;1,h); 
{h} < (s?;h,r) —h; 
endif 5 
The convex hull is {h}. The points / and r are with respectively the 
smallest and the largest abscissa. In other words they are the left-most 


and the right-most points. And k is the furthest point with respect to 1 
and r. 


Let p = {pi,p2,..-,;Pn} be a set of n generators in v-dimensional 
space, the coordinates of which are (aj), i = 1 to dand j = 1 ton. Then 
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the Delaunay tessellation in d dimensions which spans p is generated 
by, 


for j = 1 ton do 

{q} © a = (3, 0, 25)s 
endfor 
he c(q); 


project all the lower n-faces of c(q) parallel to the d‘® axis on to 
the original n-d space; 


Okabe et al (1992) give a good review of algorithms for generating 
VT’s. Given the set of generator points {p;}, i = 1 to n, a brute force 
albeit simple method generates for all i and j from 1 to n the (n — 1) 
half planes h(p;,p;), 1 <j <n, i # Jj, and then proceeds to construct all 
V(p;i) of the VT from their common intersections. 


On the other hand, the following Algorithm 9 is the quaternary in- 
cremental method whose inputs comprise the (n—3) generators p;, i = 4 
to n, where all the p; are in s = {(z,y)|0 < x,y < 1}, and three additional 
generators p, = (0.5,0.5(1 + 3V2)), ps = (0.25(2 — 3V6), 0.25(2 — 3V2)) 
and p3 = (0.25(2 + 36), 0.25(2 — 3V2)). 


Algorithm 9 Quaternary incremental method. 
k + min(k;) such that k; € St and n < 4*; 
Sij <= (i .? 1/2" ~~ 19/2"); 
construct a quaternary tree; 


scan leave buckets from left to right, top to bottom, put generators 
in buckets; 


{pi} do quaternary reordering on the generators; 
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Y ¢ construct Voronoi for p;, p2 and p3; 
for i= 4 to n do 
repeat + 1; 
while repeat do 
find pq such that d(po, pe) = min; d(p;, pe); 
if d(pm, pe) < U(pa, pe) then 
repeat < —repeat; 
elseif pm < pa do nothing; 
else 
repeat <— repeat; 
endif 
endwhile 
{wi1, wiz} < the intersections between the perpendicular 
bisector of pyp, and O(V(p;)), 1 <i<@-1; 
QO € construct the boundary of V(p;) formed by these Wyrwy7; 
{Ve} & (Ve-1} VO) — {V|V € {Ve_-1}, V is within O}; 
endfor 


(V} = Vai 
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3 T T T 


0.5; 4 


-1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 


The construction of s;; is such that s;; = (¢—1,4/2*] x (j —1,7/2*], 
where i,j € St, i and j are 1, 2, ... up to 2*. The nearest neighbour 
search finds from ¢ generators, p1,...,p¢, and {Ve_1}, the generator point 
Pm such that d(pm, pe) < d(pm, pi) for alli = 1 to £,i A mandi F £; p; are 
all generators adjacent to V(p;). The boundary growing procedure gives 
the sequence of boundaries {NO}. The additional generators mentioned 
in the procedure is graphically shown in Figure 1. 


If we draw a horizontal line through each point in two dimensions, 
or a horizontal plane each point in three dimensions, we divide the 
space into slabs the line segments within which do not intersect one 
another. In three dimensions, and for a polygonal model of porous 
media, we can for instance divide the model into such slabs, then find 
the effective cross sectional area of each slab, and then determine where 
the bottle neck to the flow occurs within the media. Whether this would 
produce the correct determination of the pressure of flow across the 
material is another matter because flows through porous media may be 
governed by the combination of the various tortuous paths through the 
pores, the interrelationship of which can be complicated. 


These slabs provide another method in finding the area, or in three 
dimensions the volume, of each cell of the tessellation. Here the cell is 
divided into slabs, and then the area or volume of each section calcu- 
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lated. 


In three dimensions, rotation around the z-axis is done by 


1 0 0 
0 cos@ —siné@ 
0 sin@ cosé@ 


around y-axis by 


cos@ OQ. siné 
0 1 0 
—sin@ 0 cosé 


and around z-axis by 


cos? —sin@é 0 
sin@ cosé@ QO 
0 0 1 


The following programme finds the area and plane parameters of a 
face from a matrix containing the list of the coordinates of the ordered 
vertices. Its synopsis is (a,p)=findfarea(v), where a is the face area, p 
the list of the plane parameters and v the coordinate matrix. 
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% findfarea.m, finds face area, (c) Kit Tiyapan, February, 2001. 
function [fca, ppr] = findfarea(odvc) 
nvthf=size(odvc,1); 
stp=floor (nvthf/3) ; 
ndpt =1+stp; 
rdpt =1+2*stp; 
nmvc= cross((odvc(ndpt, :)-odvc(1,:)), (odvc(rdpt,:)-odvc(1,:))); 
clov =[odvc; odvc(1,:)]; 
xsq =nmvc(1,1)*nmvc(1,1); 
ysq =nmvc(1,2)*nmvc(1,2); 
zsq =nmve(1,3)*nmvc(1,3); 
znmve =sqrt (xsqtysqtzsq) ; 
nmmmv =nmvc/znmvc; socp =0; 
for i=i:nvthf, 
socp =socptcross(clov(i,:),clov((i+1),:)); 
end 
fca =(abs(dot(nmmnmv, socp)))/2; 
nvct =ones(3,1); 
crdm =[odvc(1,:); 
odvc(ndpt,:); 
odvc(rdpt,:)J]; 
aprm =det([nvct, crdm(:,2), crdm(:,3)]); 
bprm =det([crdm(:,1), mvct, crdm(:,3)]); 
cprm =det([crdm(:,1), crdm(:,2), nvct]); 
dprm =det([crdm(:,1), crdm(:,2), crdm(:,3)]); 
ppr =L[aprm, bprm, cprm, dprm]; 


Real problems are algebraical or analytical whereas computer sim- 
ulations are arithmetical and numerical (cf Knuth, 1997). In between 
these two lie computer programmes. Therefore the latter are mapping 
from the analytical to the arithmetical world. There are different ways 
to solve a problem analytically, and there are different ways numeri- 
cally. Therefore numerical study by simulation is an endless pursuit. 
Several programmes given by Tiyapan (2004) do the same job, but each 
one does it differently. 


Hamiltonian’s are essentially a function that maps N bodies pair- 
wise to account for all pairs of their mutual interactions. They occur 
in most fields where there are particles interacting with each other. In 
quantum mechanics, for example, the time independent Schrédinger 
equation for a system of N particles interacting via the Coulomb inter- 
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action is HY = &Y, where the Hamiltonian is 


n> fio Bey (1) 
= \ 2m * 2 Areg|ri —75\’ 


t= 1 92 


where Y is the N-body wavefunction, z the charges of the individual 
particles and E the energy of either the ground or an excited state of 
the system. Similarly from Rushbrooke and Morgan (1961), using their 
notations, the Hamiltonian of the Ising problem is 


H= 275° 35) — gBHY- 3), (2) 
(iJ) (a) 


where £ is the Bohr magneton, g the gyromagnetic ratio and J the 
magnitude of the exchange interaction. 


When simulating such systems, the number of pairwise summa- 
tion terms can be reduced by half because they represents a symmetric 
matrix. The total number of terms is thus reduced from N(N — 1) to 
N(N — 1)/2 (cf Wray et al, 1983). Because i # j, all the diagonal com- 
ponents of the matrix are excluded, which makes the number of pairs 
n?—n =n(n—1). Ina programme, this is equivalent to two if state- 
ments, one embedded within the other in the form i(j(-)), where the 
index i runs from 1 to (n—1) and j from (4+ 1) to n; this has been 
discovered from experience, as can be seen by comparing the present 
work with the earlier one (Tiyapan, 1995, X). 

Convex Hull 


The set of extreme points EF in some superset S is the smallest 
subset of S such that the convex hulls of both £ and S are identical. 
Extreme points never lie in a triangle. 


The Graham’s scan algorithm (cf Preparata and Shamos, 1985) po- 
sitions itself in the midst of the points and then scans around in one 
direction. It determines three points in turn, and rejects a mid point 
among the three if the angle made there is reflexive, i.e. @ such that 
a > a in the anti-clockwise direction. In other words, an angle is a 
right turn if it is reflexive; it is a left turn otherwise. The algorithm is 
shown in Algorithm 10. Here {h} is a stack which contains the points 
on convex hull 


Algorithm 10 Graham's scan. 


Vaen Sryayudhya, Editor April 2006 137 


Tyabandha Journal of Arts and Science Vol. 3, No. 2 
ied; 
jo 2; 
k¢3; 
while there exist unprocessed points do 
if Zpyp;p, > 7 then 
IG 
ici-I; 
else 
tj; 
jk; 
k¢ek+]; 
endif 


endwhile o 


The algorithm starts from a known extreme point. It can also start 
from two points which are known to be extreme, in which case the 
space is divide into an upper- and a lower hulls. Similarly to the 
Graham’s scan, Jarvis’s march finds one extreme point after another as 
it wraps a line around the convex hull. The Quickhull algorithm in two 
dimensions is described in Algorithm 11. 


Figure 11 Quickhull in two dimensions. 
I + the point with the smallest abscissa; 
y ¢ the point with the largest abscissa; 


{s} < all points above Ir; 
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{s} € all points below Ir; 
while there remains an unprocessed s; in {s} do 
h © point in s; which maximises the area of Ahir; 
reject points bound by Ahir; 
{s} < all points outwards from 1h; 
{s} < all points outwards from rh; 


endfor a 


The divide-and-conquer algorithms divides a problem into sub- 
problems, finds the convex hull for each one of them and then merges 
these together by finding the convex hull of convex hulls. Algorithm 
12 finds the convex hull of two convex hulls. 


Algorithm 12 Convex hull of convex hulls 

p< one point in hi; 

if p is also in h2 then 
h © scan around from p, and merge hy and hz; 

else 
(u,v) + points on hz such that Zupu is maximised; 
ce the chain from u to v which is furthest away from p; 
h & merge cand hi; 


endif 7 


There are also dynamic algorithms for finding the convex hulls. 
In this case the input is online and one can not look ahead at the 
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input. This kind of algorithms may be useful for online applications, 
for example the traffic control in real time. Two things can happen in 
such dynamic algorithms; points are inserted or points are deleted. 


The gift-wrapping methods, which are similar to the Jarvis’s 
march, can be extended to the general d dimensions. Here a point 
is beneath a facet if it is on the same side of it as the hull, otherwise it 
is beyond it. Also, with respect to a point p, if p is beneath all facets 
that contain v then v is concave; if p is beyond the same then v is reflex; 
otherwise v is supporting. 
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